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ABSTRACT 

Numerical  methods  are  given  for  finding  unknown  functions  contained 
in  ordinary  differential  equations  when  a  solution  of  the  equation  is  known. 
These  are  iterative  methods  giving  a  best  fit  to  the  £_  norm  to  the  known 
solution,  which  may  contain  random  errors.   A  stability  requirement  on  the 
numerical  methods  is  proved. 
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1.   Introduction 

This  paper  considers  an  inverse  problem  in  differential  equations, 
that  of  identifying  ordinary  differential  equations  given  their  solution. 
Such  a  procedure  may  be  of  use  as  an  aid  in  building  mathematical  models 
to  describe  observed  phenomena. 

Since  many  different  equations  may  have  the  same  solution,  the 
problem  is  ill  posed  and  cannot  be  solved  in  this  generality.   It  is 
therefore  necessary  to  restrict  the  problem  to  cases  where  the  form  of 
the  equation  is  given  but  it  contains  arbitrary  functions.   The  goal  is 
then  to  identify  these  functions.   This  puts  the  problem  in  a  form  concrete 
enough  to  be  amenable  to  solution  by  numerical  methods.   For  example,  an 
attempt  to  verify  the  inverse  square  law  of  gravitation  might  start  with 
the  equation  (in  one  dimension) 

x"  (t)  =  f[x(t)] 

and  using  suitable  data  for  x  would  identify  the  function  f  numerically  to 
give  tabular  values  of  the  function 

f(x)  =  k/x2. 

In  order  to  simplify  the  problem,  it  is  assumed  that  each  unknown 
function  is  a  function  only  of  the  dependent  variable  or  one  of  its 
derivatives.   This  avoids  the  more  difficult  problem  of  identifying  functions 
of  two  or  more  variables. 

If  f  represents  the  unknown  functions  to  be  determined,  let 
E(f)  =  J  |xf  -  x| |  ,  where  xf  is  the  solution  to  the  equations  (the  initial 
conditions  are  ignored  for  the  moment)  and  x  the  observed  data.   Iterative 


methods  are  used  to  minimize  E  over  all  f  for  which  the  solution  x  exists. 
In  general  only  a  local  minimum  can  be  found  but  by  adding  a  "regularization" 

?  9 

term,  of  the  form   [c   f    +  c.   f    1,  to  the  error  functional  this 

o  '  '  ,.  '  '     1 '  '  _  '  '    . 

problem  is  overcome.   If  the  observed  solution  x  contains  random  errors, 

this  also  has  the  effect  of  smoothing  f. 

To  solve  the  problem  numerically,  each  f  must  be  replaced  by  a 

n 

finitie  element  approximation  f(x)~f  (x)  =  Y  q   d>.   (x)  where  the  {<b.   } 

n     .L-      i  Ti,n  Ti,n 

i=l 

are  a  given,  linearly  independent  set  (in  the  examples  given  they  are 
splines).   It  is  shown  that,  under  certain  conditions,  the  f  giving  the 
optimal  value  of  E  approaches  that  obtained  by  minimizing  over  all  f  as 
n  increases. 

The  stability  of  multi-step  methods  for  solving  differential 
equations  is  considered  from  the  point  of  view  of  estimating  not  the 
difference  between  the  computed  and  the  exact  solutions  for  a  given  equation, 
x  and  X  respectively,  but  between  equations  which  have  x  and  X  as  their 
solutions. 


2.   Solution  Methods 

Linear  problems  can  be  reduced  immediately  to  a  corresponding 

system  of  linear  algebraic  equations.   For  example,  a  system  of  n  variables 

dx 

—  =  Ax  has  the  solution 

dt 


x(t)  =  eAtx(0)  =  (eVx(O) 


=  $t  x(0)        (say) 


(1) 


if  x  =  x(i)  then 

x.+1  =  *  x..  (2) 

Taking  a  set  of  n  of  these  equations  yields  the  matrix  equation 

[x  ,  x  ,  ...,  x  ]  =  $[x  ,  x  ,  ...,  x   ],  (3) 

1   Z        n       o   i        n-l 

from  which  $,  and  hence  A,  may  be  found,  provided  the  last  matrix  is  non- 
singular.   $  can  only  be  found  to  within  a  similarity  transformation.   For 
the  case  where  the  observed  x's  certain  random  errors,  this  problem  has 
been  extensively  studied  (Lee  [6]). 

For  the  remainder  of  this  paper  only  nonlinear  equations  will  be 
discussed.   For  simplicity  only  a  single  equation  containing  one  unknown 
function  is  considered,  although  the  results  may  be  extended  to  systems. 
Usually  it  will  be  sufficient  to  consider  the  special  equation 


|f  =  f(x(t))  (4) 


The  problem  may  be  solved  by  minimizing  the  error  functional 
E(f)  =  ||xf  -  x||2  (5) 


where  x(t)  is  the  observed  solution  and 


xf(t),  o<t<T,  is  the  solution  of 


f  -f(x(t» 


(6) 


x(0)  =  x 


for  all  pairs   f,  x   for  which  a  unique  solution  is  obtained.   (x^  should 

o  n  f 

be  written,  more  precisely,  as  x,f    *  but  no  confusion  is  caused  by 

o 
dropping  the  x  ) .   The  norm  used  is  given  by 


|y||2  =  L  ^y(t)}2  dt 


(7) 


In  order  to  solve  the  problem  numerically,  it  must  be  para- 
meterized in  the  form 


n 


f(x)  =  I   q,  <fr.(x) 
1-1 


(8) 


where  {<}>.}._,  is  a  fixed  linearly  independent  set. 


n 


Let  some  set  of  n  parameters  {q. K  ,,  with  n  fixed,  together 

with  the  initial  condition  x(0)  =  x  form  the  set  over  which  otimization 

o 

is  to  take  place.   If  x  is  equated  with  q  ,  a  set  of  n  +  1  parameters 
r  o  ^o 

{q.}._  must  be  found.   The  discrete  problem  may  be  stated  as 


minimize  E(q)  =  |  |x  -  x|  | 

I  I  I  1 2   rT  2,^N  , 

where  ] |y |    =  J   y  (t)  dt 


(9) 


r 


dx 


i=l 
x(0)  =  q^ 


(10) 


V. 


If  the  function  E  is  convex  then  its  minimum  can  be  found  by  solving  the 
set  of  equations 

VE(q)  =  0. 

Lack  of  convexity  may  cause  convergence  to  a  local  minimum,  also,  if  the 
solution  is  not  unique,  difficulty  may  be  experienced  in  obtaining 
convergence  from  the  iterative  method  used.   To  avoid  these  difficulties, 
a  regularization  term  of  the  form  aW[f]  may  be  added  to  E,  where 

W[f]  =ci|f||2+c,||f'||2,  c  ,  c   >  0  (11) 

0i  i  i i      1 1  i   ii  »   Q >   i  _ 

c   and  c.  not  both  being  zero  and  a>0.   Norms  on  the  function  f  are 

°      1  k 

i  i  i  2    r        2 
given  by  i|g||   =  J   (g(x)}  dx,I  =  [a,  b] .   The  effect  of  this  is  twofold, 

a. 

-   2      n 


if  a  is  large  enough,  E  (q)  =   x  -  x    +  aW[T  q.  d>.]  is  a  convex 

a         q  j  i  i  i 

~         i=l 

functional  and  secondly  if  the  data  x  contain  random  errors,  it  will  have 
a  smoothing  effect  on  the  solution. 

To  obtain  the  derivatives  of  E  with  respect  to  q  , 

H  .  I1  [x(t) .  ;(t)1  n<«  dt  (12) 

i  i 

The  quantities  —  satisfy  a  differential  equation  obtained  from  the 
i 

original  equation  by  differentiating  it  with  respect  to  q . . 

(f(x(t)) 

(13) 

3x   9q  .    9q  t 
ni     ni 


9      /dx, 

9q.     Mtj 
l 

d 
dqt 

9f      9x 

-+9f 

Assuming  the  order  of  differentiation  can  be  changed,  and  setting 


9x  _  . 
9q~  =  6i 


and 


ff(*)    -|j     K    +j(x)    =   I   q      *'(x)  (15) 

3=1  J      J  j=l   J      J 


9f  9        n 

(x)    "  TT~     I  q,    <Mx)    "    f,(x)  (16) 


9q.  3q.   .  i  j    j  i 

i  i  j=i  * 


The   initial   condition  being   6.(0)    =   0. 

For  q    ,    the    initial    condition  at   t   =   0, 

no  ' 

dt  ^9q  ;    3x  3q     3q   »  (17) 

o         o    no 

9f 
but    /9q  =0  since  q   is  simply  the  initial  condition  x  .   The  initial 
o  ^o  o 

9x  9x 

condition  —  at  t  =  0  is  j^-   =  1.   If  ~-  e  6  ,  equations  (14)  and  (17) 

o  o  ^o  ° 


may  be   combined  as 


ft-  6i(t)  =  H  6i(t)  +  *±<*w>  (i8> 

where   d>    (x)    =   0 


with   initial   conditions 


6    (0)    =   1 
o 


6.(0)    =   0  i  =   1,    2,    ..  .,    n 

The  simplest  minimization  algorithm  is  the  "steepest  descent" 
method, 

k+1    k  ,     k 
q         -  q     +  ckP  , 

where     p   =  -  VE(q  ) 


The  scalar  c,  being  chosen  so  that 
k 


k     k 
E(q  +  cp  )  is  minimized  over  all  c>0, 


Since  E  is  assumed  convex,  at  c  =  c, 

k 

||  (qk  +  cpk)  =  0  (19) 

Since  we  are  minimizing  a  function  of  only  one  parameter,  Newton's 

method  may  be  used 

3E.   ..   CD 

(i+1)  r       (1)    3c|C 
c      =  c 


32E,   .„   (D 
c  =  c 


3c2 

3E 
To  calculate   /3c,  we  have 

T 
E(q  +  cp)  =  Jo  (xq  +  cp  -  x)2  dt  (20) 

So,  assuming  the  order  of  integration  and  differentiation  can  be  interchanged 

f  (,  +  cp)  -  2  %    (xq  +  cp  -  h   ff  dt  (21) 

3x 
The  quantity   /3c  may  be  obtained  as  the  solution  to  a  differential 

n  n 

equation,  in  the  case  where  x'  =  f(x),  f(x)  =  J  q.  cj) .  (x)  +  c  £  p.  cf).(x), 

i=l  i=l 

therefore 

d_  ,3Xv    3f  3x    3f 
dt  Sc     3x  3c    3c 

(22) 

n 

3f 


and       "3c"  ml     Pi  *i(x) 


i=l 
3x,„     ......  .      3x 


n 


Thus  "/3c  satisfies  the  same  equation  as   /gq.  except  ^  p   <j> .  (x)  replaces 

1  •   -i    X    J- 

3x.  1= 

tj).(x).   The  initial  condition  is  —        =  p  . 
l  3c   t  =  o    o 


For  the  second  derivative, 

2  T 

3  E  /      n    of   /8x\2  J4. 

-2   (,  +  cP)  -  2  /o  (_,  at  (h) 

T         A    2 
+  2  /   (x  ,    -  x)  ^-f-  dt 

oc 

a2x   2 

An  equation  for    /3c  may  be  found  or  the  assumption  made  that 

x  ,    is  close  to  x  and  the  second  term  ignored.   Thus  only  one  or  two 
q+cp 

differential  equations  need  be  solved  at  each  iteration  of  the  calculation 
of  c. 

Other  algorithms  may  be  obtained  using  different  choices  for 
p  .   For  instance  Newton's  method 

pk  =  -  {Hk}_1  VE(qk)  (24) 


where  H  is  a  matrix  with  components 

H   -£ 

Id    3q±  9qj 

■ 2  /  ¥■  ¥  dt  <25> 

■'o  3q  .  3q. 
i   J 

T      A    2 

+  2  f   (x  -  x)  ^     dt. 
Jo  8qi  3q 

32x 
The   equation   for  /3q.    3q . ,    for   the  case  x'   =   f(x)    is 

i        J 

d      x2  32f    x      x 

— —   o     .    =   — —   0  .     0  . 

dt      U        3x2      1     J 

3x      ij 


(26) 


+  f-   *.(x)    6. 


+  fc*j(x)  {i 


,     .   _  3x   s2  9  x 

where  o.  -   —  ,  6..  =  7 — . 

1    dq±'      ij    3q±3q;j 

2 
This  however  requires  the  solution  of  about  n  /2  equations,  if  x  -  x  can 

be  assumed  small  the  second  term  may  be  neglected.   This  latter  method 

(the  Newton-Gauss  method),  works  well  in  cases  where  x  contains  very  little 

random  error  and  a  closely  approximating  x  can  be  found. 

A  comparison  of  gradient  methods  for  parameter  estimation  problems 

is  given  by  Bard  [2] . 


10 


3.   Stability  Analysis 

Consider  our  test  equation  dx/dt  =  f(x)  and  suppose  some 

numerical  solution  {X  },  i  =  0,  1,  2,  . . . ,  n  is  obtained  where  X  is  the 

1  i 

solution  at  t  =  ih.   Let  X(t)  be  any  function  taking  the  values  X  at 

i 

t  =  ih  and  suppose  F(X)  is  some  function  satisfying  dX/dt  =  F(X).   Let 

x.  =  x(ih). 

l 

The  usual  stability  analysis  is  concerned  with  the  behaviour  of 
*L  -   xi  for  increasing  i  (here  called  "stability  in  x") .   Since  we 

are  concerned  with  the  inverse  problem  it  is  necessary  to  find 

numerical  methods  for  solving  equations  for  which  it  can  be  shown  that  F 

is  close  to  f.   As  the  functions  are  evaluated  only  at  the  points  {x.},  the 

behaviour  of  F(x.)  -  f(x.)  is  analyzed  for  increasing  i.   Roughly  speaking, 

a  method  may  be  said  to  be  "stable  in  f"  if  F(x.)  -  f(x.)  remains  bounded. 

Our  analysis  is  restricted  to  multistep  methods. 

The  general  k-step  method  of  order  q  may  be  expressed  in  the  form 

k 

T   [a.  x.  .  +  he.  f(x.  .)]  =  S.  (27) 

where     S.  =  C  hq+1  x.(q+1)  +  0(hq+2) 
i  l 


and 


k    (    n\1  (    \\^~x 

siQ     q!   3      (q-D!   j 


k 
Let       p(C)  =   I     a   r  J  (28) 

j=o 

k 
a(0    =   I     6  C   J  (29) 

j=o 

The  condition  that  the  method  be  weakly  stable  in  x  is  that  the  roots  of 
the  equation 


P(0  -  0  (30) 


11 

lie  within  the  unit  circle  or  are  on  the  unit  circle  (see  Gear  [4]  §8.3. 

(The  method  is  strongly  stable  (in  x)  if  the  roots  lie  within  the  unit 

circle  except  for  the  root  £  =  1) .   We  show  that  a  similar  condition  on  the 

roots  of  a(£)  =  0  is  required  to  ensure  stability  on  f. 

Assuming  that  equation  (27)  is  used  as  a  corrector  and  that 

sufficient  iterations  of  the  corrector  are  made  then  this  equation  holds 

for  the  computed  X  except  that  the  truncation  term  is  replaced  by  zero 

and  some  roundoff  error  is  incurred. 

k 

y   [a.  X.  .  +  h  3.  f(X.  .)]  =  R.  (31) 

j=o   J   1_J      J     1_J      x 

Equation  (27)  also  holds  for  the  solution  X  of  X'  =  F(X). 
k 


y   [a.  X.  .  +  h  3.  F(X.  .)]  =  T\ 


(32) 


J=° 


where     T.  =(hq+1  X.  (q+1)  +  0(hq+2). 


If  we  define  AF.  =  AF(X.)  =  F(X.)  -  f(X.),  subtracting  equation  (31)  from 

equation  (32)  gives 

k 

I     h  6.  AF.  .  =  T.  -  R.  (33) 

3  =  0     l-l-i 

This  equation  determines  the  difference  between  F  and  f  at  the 

solution  points  {X.}.   The  analysis  must  be  restricted  to  the  set 

{AF(X.)}.   A  method  will  therefore  be  said  to  be  stable  in  f  if  a 
i 

perturbation  AF  will  produce  subsequent  changes  which  do  not  increase 
step  to  step. 

The  general  solution  of  the  difference  equation  (33)  may  be 
found  by  taking  any  particular  solution  and  adding  to  it  any  linear 
combination  of  solutions  of  the  corresponding  homogeneous  equation 


12 


)      h   B.    AF.     .   =   0.  (34) 

A  particular   solution  of    (33)    is 
i      (T.    -  Rj 


AF.  =  7  — i— - — ^—  y. 


where  y.  satisfies  (Henrici  [5]  §5.2-1) 


i 


y.  =  0  o<i<k-l, 

J  l  

X    h  3.  y.  .  =  [ h  for  i=k 

j=o    3      X  J    So  for  i>k 

The  general  solution  of  equation  (33)  is  therefore 


(35) 


k  i   (T.  -  R.) 


af.  =  y  c.  z.  .  +  y  — i 
-,=1  j     »j   -|=k 


J-y,  4,  (36) 


h     'i-j 


where  the  sequences  {z.  . }  (j  =  1,  2,  . . . ,  k)  are  the  set  of  linearly 

i»  3 

independent  solutions  of  the  homogeneous  equations,  given  by 


z.  .  =  r.  ,  (37) 

i.J    3 

where  r.  is  a  root  of  the  equation 
J 

k        . 
a(r)  -J   6.  r  J  =  0,  (38) 

j=o   J 

or,  in  the  case  where  a(r)  =  0  has  a  root  of  multiplicity  m  then 

i   .  i   . 2  i       . m—  1  i   ,,,  .     ..      _ ,  ,    r- ,  , 

r  ,  lr  ,  i  r  ,  ...,  i   r  wxll  be  solutions.   Since  y.  satisfies  the 

homogeneous  difference  equation  for  i>k,  y   is  a  linear  combination  of 

the  z.  .  for  i>k.   This  proves  the  following  result. 


13 

Theorem 

A  multi-step  method   is   stable   in  f   if  and  only   if   the  roots  of 
the   equation 

o-U)    =   0 

lie  within  the  unit  circle  or  are  simple  on  the  unit  circle. 

A  method  may  be  said  to  be  strongly  stable  in  f  if  all  the  roots 
of  o(0   =  0  lie  within  the  unit  circle  and  weakly  stable  if  any  simple 
root  lies  on  it. 

The  well  known  result  of  Dalquist  [3]  that  a  k-step  method 
which  is  strongly  stable  (in  x)  is  of  order  at  most  k+1  (for  k  even  a 
weakly  stable  method  of  order  k+2  is  possible)  may  be  extended  to  cover 
the  case  where  stability  in  f  is  also  required. 

Theorem 

A  k-step  method  which   is   stable  both   in  x  and   in   f  has   order 
at  most 

k  if  k  is   even 
k+1   if  k  is  odd 

Further   if   the  method   is   of  order  k+1,    it   is  weakly  stable   in  both  x  and   f , 

Proof 

The  proof  is  similar  to  that  used  to  obtain  Dalquist 's  original 
result.   The  condition  that  the  method  be  of  order  r  is 

f^-(0  +  o(0   =  0  (U-l)r).  (39) 


14 


Making   the  transformation 

r   =  2±£         z    =   ill 
*       1-z   »  5+1 

The  interior  of   the   unit    circle   in   the   £-plane  maps   into   the  left 
half  of   the  z-plane 


Let  R(z)   +    (^)k  p(^f) 


and  S(z)   =   (^)k  a(±±§-) 

z  1-z 


(40) 


These  are  both  polynomials  of  degree  k  in  z.   Except  at  z  =  1  and  E,  =   -1 
the  roots  of  R  and  p,  and  of  S  and  a   correspond.   Since  p  and  a  are  stable, 
£  such  that  z  =  1  is  not  a  root  of  p(0  =  0  or  o(0  =  0.   The  roots  corres- 
ponding to  E,   =  -1  in  the  z-plane  are  at  infinity  and  so  decrease  the  order 
of  R  or  S. 

Equation  (39)  maps  into 

+  S(z)  =  0(zr).  (41) 


log  {(l+z)/(l-z)} 

Following  Gear    [4]      10.2,    let  R(z)    =   a     +  a   z  +   ...   +  a  z    . 
Since    E,  =   1    is    a   simple   root   of  p(£)    =0,    z   =   0    is    a   simple   root   of  R(z)    =    0. 
Therefore,    a     =   0  and  a     ^   0. 

Since  p,    and  hence  R,    is    a  real   polynomial,    the  roots  of  R(z)    =   0 
are   either   real,    x   ,    or  occur   in  conjugate  pairs,   x      ±   iy     with  non-positive 
real  parts.      Hence 


R(z)    =   air      (z-x    )tt      (z   -   x      -   iy   )    (z    -   x     +   iy   ) 

\i  yv  vv  v  v 

o  9  2 

=   a7T      (z+|x)tt      (z     +2z|x|+x'''y      ). 

y  '    y       V  v '  v        v 
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Therefore   the  signs  of   the  non-zero   a.   are  all   the  same.      Without   loss  of 
generality,    take  a   >0,    then  a.>0,    i>2. 
Similarly   if 

S(z)   =b     +  b-    z  +'...+  b.z   ,  (42) 

o  1  k 

then   all   the  b.   have  the  same  sign  since   the  roots  of   S(z)    =   0   also   have 
non-positive  real  parts. 

Tf  ?. y    c     z2^ 

Lt  log  {(l+z)/(l-z)}   y£Q  °2y 

it  may  be  shown  that  (Gear  [4]  §10.2) 


(43) 


Co  =  \   '  C2y  K   °'  y^- 

Then    Ml) =  y  a  z™-1  y  c   z2y 

Then      log  {(l+z)/(l-z)}   I     am  Z    *  °2y  Z   ' 

m=l        y=o 

00 

=  I      rn  z  ,  say 
n=o 


(44) 


In  order  for  the  method  to  be  of  order  p, 

00  If 

y  r  zn  +  I   b  zn  =  0(zP),  (45) 

^   n      L     n 
n=o        n=o 

so  that   r  =  -b   for  n  =  0,  1,  2,  ....  p-1 
n     n  '  '  r 

and  hence  the  coefficients  r  must  all  have  the  same  sign  for  0  up  to  p-1. 
Equating  coefficients  in  equation  (44) 
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r     =   c     an 
o  ol 


r     =  c     a0 
1  o      2 


r0  =   c      a„  +  c„   a, 
2  o      3  2     1 


r.=   c      a,    +  c0   a0 
3  o      4  2     2 


rk 


c0  a.         +  c.    a,     _+...+  c1         a.   +  c.  a.  k  even 

2      k-1  4      k-3  k-2      3  k  1 


c      a,     ..    +  c.    a.     _+...+  c,     „   a,    +   c,     .  a„      k  odd 
2      k-1  4      k-3  k-3      4  k-1    2 


Since  c      and  a,    are  both   >  0,    r     >   0.      Therefore  r, >0   for  r  =   1,    2,    ....    p-1 
o  1  o  1—  r 

However,    since  a,    0,    a.>0   for   i>2   and   c„   <0  for   y>l,    r   <0   for  k  even  and 
1  l—  —  2y  —         k 

r  < 0   for  k  odd.      Therefore   for  k   even,    the  method  has  maximum  order  k, 

K, 

and  for  k  odd  the  method  is  limited  to  order  k+1  by  the  requirement  for 
stability  in  x  . 

When  the  method  is  of  order  k+1,  (and  k  is  odd),  r  must  be 

K. 

zero  and  so  a0  =  a,  . . .  =  a,  n  =  0.   Since  a  is  also  zero,  R(z)  is  an 

Z      H  R— 1  O 

odd  polynomial.   Further  r..  =  r_  =  ...  =  r,  =  0,  so  b.  =  b„  =  ...  =  b  =  0 
and  consequently  S(z)  is  an  even  polynomial.   Therefore  if  x  +  iy  is  a 
root  of  either  R(z)  =  0  or  S(z)  =  0  then  so  is  -x-iy,  but  the  roots  of 
these  equations  must  have  non-positive  real  parts  to  satisfy  the  stability 
requirements  and  so  all  the  roots  lie  on  the  imaginary  axis.   The  root 
of  both  p(0  =  0  and  o(0    =   0  therefore  all  lie  on  the  unit  circle. 
Example:   If  B  =  1,  3  =  $     =  3  =  B,  =  0,  then  all  the  roots  of  o(0  =  0 
are  zero,  and  (directly  from  equation  (33)), 

Ti  -Ri 
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This  yields   the  fourth  order  four-step  method, 


j 

> 

0 

1 

2 

3 

4 

J 

-25/12 

4 

-3 

4/3 

-1/4 

3 
j 

1 

0 

0 

0 

0 

for  which   the   roots  of  p(£)    =   0  are  1,    0.3815,    0.2693    ±  0.04920i. 
method   is    therefore  strongly   stable  both   in  x  and   in   f. 


The 
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4.   Test  Problems 

In  the  first  example,  an  additional  test  was  made  in  which  the 
observed  data  were  perturbed  by  adding  random  "noise".   This  was  mainly 
to  test  the  effect  on  the  convergence  rate. 

Example  1 

This  is  the  simplest  possible  example,  taking 

dx    cr    \ 

dT  =  f  (x) 

and  input  data  of  the  form 
-1 


x(t)  = 


t+2 


specified  on  the  interval  [0>10]  at  discrete  points  distance  .05  apart 
This  will  give  the  solution 


f(x)  =  x2 

on  the  interval  (determined  from  the  range  of  the  observed  data)  [-   /2,    -   /12] 
A  smoothing  factor 

i  i  i  i2       ii   i  i2 

a [c    : |f  1 1    +  c1  : |f '  1 1  ] 

-2 
was  introduced  with  c  =  0,  c  =  1  and  a  initially  set  to  10   and  reduced 

-9 
by  a  factor  of  10  at  each  iteration  until  a  minimum  of  10   was  reached. 

The  function  f  was  approximated  on  the  interval  [-  /2,  0]  by  a  natural 

spline  with  8  equally  spaced  knots.   (Natural  splines  are  discussed  in 

Ahlberg,  Nilson  and  Walsh  [1]). 

At  each  iteration  the  value  of  the  error  functional  and  the 
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%     norm  of  its  gradient  were  produced.   The  programs  were  rerun  after 
applying  a  random  noise  with  zero  mean  and  a  variance  of  .01  to  the  input 
data.   It  may  be  noted  that  in  all  cases  the  function  was  better  identified 
in  the  middle  of  the  range  where  more  data  points  were  avilable  than  on 
the  edges . 

Table  1.   x'  =  f(x)  -  no  noise  added  to  the  observed  x. 


Iteration 

|  | VE  |  | 

E 

1 

6.77E-00 

9.99E  02 

2 

5.56E-00 

5.14E-02 

3 

4.91E-00 

2.60E-02 

4 

1.04E-01 

1.22E-04 

5 

7.09E-02 

8.09E-06 

6 

2.54E-03 

6.00E-08 

7 

7.76E-05 

2.38E-09 

8 

6.89E-06 

3.33E-11 

X 

-0.50 
-0.45 
-0.40 
-0.35 
-0.30 

-0.25 
-0.20 
-0.15 
-0.10 
-0.05 


-0.00  0.00017 


computed 

exac 

0.24993 

0.25 

0.20250 

0.2015 

0.16000 

0.16 

0.12250 

0.1225 

0.09000 

0.09 

0.06250 

0.0625 

0.04000 

0.04 

0.02250 

0.0225 

0.01000 

0.01 

0.00252 

0.0025 
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Table 

2.   x' 

=  f(x)  -  with 

noise  adde 

Iteration 

|ve| 

1 

6.86E-00 

1.01E-01 

2 

5.53E-00 

5.87E-02 

3 

4.29E-00 

2.72E-02 

4 

2.11E-00 

5.18E-03 

5 

1.64E-01 

1.21E-03 

6 

2.41E-01 

9.90E-04 

7 

4.32E-02 

9.38E-04 

8 

6.52E-03 

9.32E-04 

9 

1.69E-03 

9.32E-04 

X 


f(x) 


computed 

ex 

0.50 

3.93796E-01 

0.25 

0.45 

1.89431E-01 

0.2015 

0.40 

1.84711E-01 

0.16 

0.35 

1.03462E-01 

0.1225 

0.30 

8.22265E-02 

0.09 

0.25 

7.12669E-02 

0.0625 

0.20 

4.00412E-02 

0.04 

0.15 

2.53864E-02 

0.0225 

0.10 

1.01355E-02 

0.01 

0.05 

-5.73504E-02 

0.0025 

0.00 

-8.29363E-01 

0 

Example  2 

This  example  uses  a  system  of  two  equations  with  two  unknown 
functions.   Starting  with  the  vector  equation 


d2x 


dt' 


=  f(x), 


and  splitting  it  into  horizontal  and  vertical  components  x1  and  x~, 
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we  try  the  equations 


d  xl  2     2 

-I"  =  fl  (xl}  f2  (X1  +  X2  ) 
at 


d  x?  2     2 

-T  -   fl  (X2}  f2  (X1  +  X2  > 
at 


If  the  original  equation  represents  an  inverse  square  law  then  we  should 
obtain  for  the  functions  f..  and  f 

f^x)  =  -kx 

f    (   n    1   "3/2 
f2(x)  =  ^  x 


where  k  is   an  arbitrary   constant. 

The  non-uniqueness   of   f     and   f„   is  handled  by   introducing   a 
smoothing  term 

atc  I         I  lf  •  I  I2  +  ci         I         I  lf  *l  |] 

lq  L  ll-.ll  L  II  ||J 

i=l,2  X  L   i-1,2 

into    the   error   functional,   with   c     =   c     =1.      Again  a  was   allowed  to  vary 

o    1 

-2  -9 

from  10   down  to  10   when  noise  was  not  added  to  the  data  and  down  to 

-  fi 
10   when  it  was.   The  exact  solution  is  scaled  so  that  the  mean  value 

for  f„(x)  is  the  same  as  the  computed  one. 


22 


Table   3. 

Iteration 

1 

5.72E2 

7.19E-1 

2 

5.66E-1 

1.70E-2 

3 

1.96E-1 

9.77E-4 

4 

1.24E-2 

1.91E-4 

5 

3.29E-2 

1.3  7E-4 

6 

3.73E-2 

1.01E-4 

7 

3.74E-2 

6.66E-5 

8 

3.29E-2 

4.00E-5 

9 

2.33E-4 

1.02E-7 

fx(x) 
computed  exact 


f2(x) 
computed  exact 


-0.23 

-0.06 

0.12 

0.29 

0.47 

0.65 
0.82 
1.00 
1.17 
1.35 


0.1375 

0.0457 

-0.0965 

-0.2453 

-0.3923 

-0.5396 
-0.6868 
-0.8340 
-0.9837 
-1.1113 


0.1924 

0.0477 

-0.0971 

-0.2418 

-0.3865 

-0.5313 
-0.6760 
-0.8207 
-0.9655 
-1.1102 


1.00 

1.2006 

1.2153 

1.14 

0.9870 

1.0098 

1.28 

0.8302 

0.8430 

1.41 

0.7111 

0.7226 

1.55 

0.6170 

0.6284 

1.69 

0.5426 

0.5530 

1.83 

0.4848 

0.4407 

1.97 

0.4414 

0.4407 

2.10 

0.4097 

0.3981 

2.24 

0.3859 

0.3619 

1.53 


-1.1137 


-1.2549 


2.38 


0.3763 


0.3309 
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